Granular systems in the Coulomb blockade regime 
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Disordered granular systems, at temperatures where charging effects are important, are studied, 
by means of an effective medium approximation. The intragrain charging energy leads to insulating 
behavior at low temperatures, with a well defined Coulomb gap. Non equilibrium effects can give 
rise to a zero temperature transition between a metallic, gapless phase, and an insulating phase. 



I. INTRODUCTION 

Charging effects modify the transaoHt properties of 
granular metals at low temper aturesautl. In magnetic 
systems, the low temperature rn^ngtatpsistance can also 
be modified by charging effectsQiJ'BLrffl. 

A successful scheme used to analyze Coulomb blockade 
in granular systems is to reduce the system to the study 
of an effective single junction, as in the pioneering study 
done in Ref. This assumes that each junction has 
an activated conductance, which is approximately cor- 
rect when the high temperature conductance of the junc- 
tion g, is low in units of e 2 /h. A single junction in the 
Coulomb blockade regime has been extensively studied 
in recent yearstil. Modern techniques allow us to study 
the behavior of a single junction for arbitrary values of 
the conductance and offset charges (see below). 

The current approach describes a single junction in 
terms of two collective variables, the charge, Q, and 
its conjugate phase variable, cA_JXhe electronic degrees 
of freedom are integrated outE3Ej. This procedure is 
motivated by a similar approach used in the study of 
small Josephson junctionalJ. The present tools allow us 
to study junctions with high values of the conductance, 
g 3> e 2 /h, and to consider the effect of gate voltages, 
and of degeneracies between different charge states. The 
model can be extended to study non e quili brium pro- 
cesses associated to the tunneling eventsOEJ, which can 
be considered as the leading corrections in an-eitHansion 



in <? _1 , where g is the intragrain conductancdiZrlia. 

In the following, we will extend these standard meth- 
ods used in the study of Coulomb blockade in single junc- 



tions to the analysis of granular materials. In order to 
do so, we will use an effective medium theory. This ap- 
proach can be viewed as an extension to mesoscopic sys- 
tems of the Dynamical Mean Field Theory used in the 
study of-ttie microscopic properties of strongly correlated 
systemsta. The model and method of calculation will be 
presented in the next section. Section III discusses results 
obtained using analytical large-N techniques-. sthich de- 
scribe well the properties of single junctionsEJEa. Then, 
we show numerical results obtained using Monte Carlo 
(MC) methods previously applied to the study of single 
junctionsc 2 ! Section IV contains the main conclusions of 
our work. 

II. THE MODEL 



A. Effective medium theory. 



We will consider a system made of metallic grains. We 
assume that the separation of the electronic states within 
each grain, Si, is much smaller than the charging ener- 
gies, E l c , or the temperature, k^T. The grains are con- 
nected by a large number of channels with small trans- 
mission coefficients, so that a perturbative treatment of 
the coupling induced by each individual channel is pos- 
sible. Then, thje-fticptronic degrees of freedom can be 
integrated outE^lialiJEa, and the system is described in 
terms of the charge at each grain, Qi, and its conjugated 
phase variable, 4>i, where [Qi,(f>i] — ie, and e is the unit 
of charge. The effective action which describes the charge 
effects in the system is, at zero temperature: 
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l <3 



(1) 



where K(t) — [7rTcsc(7rTr)] 2 (we henceforth set ks = 1)- 
We further assume that the capacitance matrix is dom- 
inated by its diagonal entries, and identify the charging 
energy of the i th grain as E\ = e 2 /2C^. The coupling be- 



tween grains i and j is onj = 2\tij\ 2 pipj, where Uj is the 
hopping matrix element between the grains (assumed in- 
dependent of momenta), and pij are the respective den- 
sities of states. One can also write — Rk / Rij , 



2 



where Rk — h/e 2 — 25, 813 SI is the quantum of resis- 
tance and Rij is the high-temperature value of the inter- 
grain resistance. 

The above model possesses a global 0(2) rotational 
symmetry, which may be generalized to an S\J(N) sym- 



metry as follows. On each site, we replace the unimodular 
complex scalar Zj = e l ^ i with an A^-component unimod- 
ular complex vector with components z,^ (// = 1, . . . , N). 
The action for the SU(iV) model is then 



P p p 

S it j dT ( d r^)(drz^) + aij JdrjdT'K(T - r') |l - z.* I (r)z w (r)z w (r')^(r') 

l<3 ' 



(2) 



We now assume that each grain is coupled to a large 
number, M, of other grains, and define the mean field 
correlation function Gi(r — r') as 

— Gi{T-T') 5^ = ^c*ij z iM (r)z* 1/ (r')) , (3) 

3 

with a = Af otij) the average coupling of a grain 

to its neighbors. Note that G(0) = 1. This allows us to 
write the action as a sum over contributions from indi- 
vidual grains, with the effective self-interaction term 

P P 

= a JdrJdr'KiT-T') [l-G^r -t^z^z^t')) 



(4) 

with a = Ma/N. Thus, the properties of the system 
can be reduced to those of a single junction embedded 
in an effective environment whose properties are to be 
determined self consistently. 



electronic levels within the grain. This effect can be ana- 
lyzed using similar techniques to those used in relation to 
the sudden ejection of a core electron in a metalcJ, and 
it leads to a modification of the kernel in the term which 
describes the tunneling in the action in Eq.(^). The de- 
cay of the interaction K(t — r') at long times becomes 
\t — t'| 2_£ , where e is a dimensionless parameter, which 
can be expressed in terms of the phaseshifts induced at 
the Fermi level by the sudden switching of the external 
potential. A similar effect is expected if the tunneling 
takes place between quasi one-dimensional systemstJ. 

In fact, for each pair of grains one should replace the 
kernel K (r) with 

K«(t) = {KT [Eir krcsc(Trrr)] 2 — * . (5) 

For the sake of simplicity, we will assume 

K{t) = ST [ttT csc^Tt)] 2 ^ (6) 

for all pairs, where E c is the average charging energy. 



B. Offset charges and non equilibrium effects. 

We can extend the model discussed so far to situations 
when there are offset charges, (ft, at each grain. These off- 
set charges lead to additional phases in the contributions 
of paths of 4>(t) with non zero winding numbers. A path 
of winding number m acquires a phase Xi = 27rm<ft . This 
phase has to be added to the effective action, Eq.(Q). 

The same scheme can also be used to include the ef- 
fect of nop equilibrium processes associated to intergrain 
tunnelingliltj. The electron which tunnels induces an in- 
homogeneous electrostatic potential, which modifies the 



C. Monte Carlo procedure 

In our effective-medium scheme, the grand partition 
function for a generic grain i with a given offset charge 
qi can be written in terms of the phase 4>i , ascia: 

oo „ 

z (Qi) = ^2 ex P(%) / exp(-5eff[^i]) , (7) 

m— — oo 

where the paths 4>i{t) satisfy in sector m the boundary 
condition 4>i(j3) = 0i(O) + 27rm. The effective action takes 
the form: 



S cS [^} = ^J q dr {^) + drj dT'K(r-T') {1 - c 0S #(r) - </>(r')] G(r-r')} (8) 

I 

with the kernnel Kir) given in Eq. (^|) It describes low- the unscreened average charging energy, E c . Here, the 
energy processes below an upper cutoff on the order of 
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correlation function G is given by 

G(t - t') = (co#(t) - 4>{t')]) 



(9) 



where the average value is calculated with the partition 
function in Eq. (|7]). Note that the action in Eq. (0) is 
similar to that corresponding to a single tunnel junction 
with dimensionless conductance a and charging energy 
E c , the only difference being the presence of the correla- 
tor G in the tunneling part of the effective action of the 
granular system. 

Different sectors (winding numbers m) in the partition 
function Z(qi) have been sampled by path-integral Monte 
Carlo for temperatures down to T = E c /50. MC simu- 
lations have been carried out by the usual discretization 
of quantum, paths into N (Trotter number) imaginary- 
time slicesc3. To keep roughly the same precision in the 
calculated quantities, as the temperature is reduced, the 
number of time slices has to be increased as 1 /T. We have 
taken TV = Af3E Cl and thus, we employed an imaginary- 
time slice At — (3/N = 1/(4_E C ). Details on this Jiiud 
of MC simulations were given in earlier publicational3'E3 
and will not be repeated here. 

^From a computational point of view, the main dif- 
ference between-the present procedure and that used in 
earlier worksE3E3 is the requirement of self-consistency 
for the correlator G(r — r'), as expressed by Eqs. (g) 
and (Q). To obtain this self-consistency, we employed an 
iterative procedure that led to convergence of the cor- 
relator in the region of parameters studied here. For a 
given set of parameters (T, a, e, and qi), we begin the 
simulations by assuming a correlator Gq{t — t') = 1. The 
correlator G\ obtained by introducing Go in the effective 
action Eq.(|J) is then employed in a new iteration to ob- 
tain G2, and so on. In our simulations, each iteration 
consisted of 2 x fO 4 MC steps (in a MC step all path 
coordinates and winding number are updated). By using 
this procedure, the correlator converges after a few iter- 
ations. In general, the number of iterations required for 
convergence increases as the temperature is lowered, but 
in fact 10 iterations were enough to find convergence in 
all cases presented here. We have checked that employ- 
ing different starting correlation functions Go(r — r') in 
our iterative procedure does not change the results. 

The conductance g{T) has been obtained from Xhc 
(self-consistent) correlator G by using the expression!!!!: 



g(T) = go 



(10) 



where go is the high-temperature conductance. 



III. RESULTS 



A. Large-iV Theory 



We now address the self-consistent single grain theory 
within the large-iV approximation, where N is the num- 
ber of components of the field z^. The unimodularity 



constraint 



N 

^=1 



1 can be imposed with a real 



Lagrange multiplier field, A(r): 







1 



(11) 



s = J (t^I^I 2 + a(M 2 -i) 




+a JdrJdr'KiT -t') {l - G(r - t^z^z^t')} 


where G(t — r') = (z* (T)z M (r')) is the self-consistent 
field due to neighboring sites. The partition function 
is then obtained by integrating exp(— S) freely over the 
complex fields {z M (r), z m (t), A(t)}. In the N — > 00 limit, 
the saddle point of the functional integral is dominated 
by configurations in which A(r) is constant. The action, 
up to a harmless constant, is 



\zp,(io n )\ 2 - A 



(12) 

where we now express all energies in units of E c . We 
define L = f3E c ; the bosonic Matsubara frequencies are 
u) n = 2irn/L. The function Q(oj n ) is the Fourier trans- 
form of Q{s) = K(s)G(s) 1 where s = E c t is the dimen- 
sionless imaginary time, and 



K(s) = 



7T 

— CSC 



7T|£| 

L 



(13) 



Thus, 



Q(0) - Q{io n ) = ds(l- cos(w„s)) K(s)G(s) . (14) 



(15) 



The correlation function G(s) is given by 



4 



wa+A + fi[Q(0)-Q(w n )]- 



and the fixed length constraint is G(0) = 1. 

We iterate the mean field equations at finite temper- 
ature to obtain A(<5, L) in the following manner. First, 
given an initial guess for the set {Q(co n )}, we find A such 
that G(0) = 1. We then recompute Q(s) = K(s)G(s) 
and a new set of Fourier components {Q'(uJ n )}. Iterat- 
ing until self-consistency, we obtain the results shown in 
Fig. |l|. The self-consistent solutions reveal a discontin- 
uous jump in A at a c , i.e. the transition apparently is 
first order. 

However, it is hardly clear that our iteration scheme 
should converge to the proper (i.e. lowest free energy) 
self-consistent solution. If the transition is second order, 
we can make progress by assuming, at criticality (L = 00, 
a = <S C ), that 



00 



cos(ws) 



uj 2 + & c [Q(Q) - Q(u)] 



(16) 
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with Q(s) = G(s)/\s\ 2 - e . Since G(0) = 1, we have that 
Q(s) ~ |s| e ~ 2 as s — > 0. For s large, we assume Q(s) ~ 
Ajsl* -2 and solve for A and S by ignoring the first term 
in the denominator of Eq. (|l6|). We find 



A = 



( N 



\47ra. 



(2 - e)ctn(7re/4) 



1/2 



(17) 



and S = he. Accordingly, we adopt a trial solution for 
Q(s) which interpolates between the small and large s 
limits: 



Q00 



e -n\s\ - e -v\s\) 



(18) 



with rj a variational parameter. 

Given this trial form for Q(s), we solve G(0) = 1 to 
obtain a c (e, 77). We then check for self-consistency, com- 



paring G(s) obtained from Eq. (16) with Q(s) \s\ c . We 



find (i) that a c is insensitive to rj for e < 0.25, but increas- 
ingly sensitive for larger values, and (ii) that there is an 
optimal value for 77, in order that the two definitions of 
G(s) agree best (see Fig. ||). 

Returning to the first order behavior found in the it- 
erated self-consistent solution, we find that this behavior 
persists for e>0.30. Furthermore, it is hysteretic. The 
solution shown in Fig. [I] was obtained by sweeping a 
up from a = 0. When a self-consistent solution was ob- 
tained, the set {Q{u>n)} was chosen as the initial con- 
ditions for iteration for the next value of a. When we 
sweep both up and down, however, we find the solution 
is hysteretic, as shown in Fig. |j| In Fig. |] we com- 
pare a c obtained assuming a second order transition, as 
discussed above, with the upper and lower critical val- 
ues obtained by iteration of the self-consistent equations 
(with L = 128). 



B. Monte Carlo simulations 

We now turn to the results of our path-integral Monte 
Carlo simulations for the conductance of the granular 
system. 

In Fig. ^ we present the temperature dependence of 
the conductance g(T), as derived from our self-consistent 
correlation function G (squares) for a — 0.1 and e = 0, 
corresponding to a Coulomb-blockade regime. For com- 
parison, wc also show MC results found for the same 
set of parameters for a single junction (circles), for which 
one finds at low T an algebraic dependence of the conduc- 
tance in the form g(T) ~ T 2 , as expected for co-tunneling 
processes. These processes do not show up in the granu- 
lar system, and in this case g(T) decreases exponentially 
as the temperature is lowered. 

Our Monte Carlo procedure allows us to study also the 
dependence of the conductance of the granular system as 
a function of the offset charge wich appears in the 
partion function through the phase Xi- This dependence 



is shown in Fig. ^ for the same a value as in Fig. || at 
a temperature T — E c /10. One finds a maximum of the 
conductance at the degeneracy point (qi — 0.5), similar 
to the well-known case of a single junction in the pres- 
ence of Coulomb blockade (solid line in Fig. |^). Again, 
we observe a decrease in the conductance for the gran- 
ular system, as compared with the single junction with 
the same effective parameters (<S,e), in the absence of 
charge-offset (qi = 0). However, close to the degeneracy 
point, we find similar values for both systems, that in 
fact coincide for qi — 0.5 within the error bars of our 
numerical results. 

The conductor-to-insulator phase transition at low 
temperatures in the parameter space (e, a) can be also 
analyzed fron MC simulations. With this purpose, we 
have studied the dependence of the conductance on a 
and T for several values of e. For a given e one expects a 
conducting regime for a larger than some critical value a c 
(see Fig. above). In our numerical analysis this means 
that we expect a crossover from a region in which g in- 
creases (large a) to another regime in which g decreases 
(small a) as T is lowered. This is shown in Fig. [?], where 
we have displayed the dependence of the conductance g 
on a for e = 0.5 at several temperatures. Each kind of 
symbols represents MC results at a given temeprature, 
with T increasing from top to bottom in the right part 
of the Figure. From the crossing point of these curves we 
can estimate a critical value a c — 0.46 ± 0.02 for e = 0.5. 

In a similar way we have obtained a c for other e values, 
and the results are presented in Fig. ||. Symbols are data 
points derived from MC simulations, and the dashed line 
is a guide to the eye. 



IV. CONCLUSIONS 

We have presented a method to analyze the conduc- 
tance of metallic granular systems in the range of param- 
eters where Coulomb blockade effects are important. The 
method is a natural extension of the techniques widely 
used to study the competition between charging effects 
and tunnel conductance in single junctions. The method 
allows us to include disorder, and becomes exact when 
the average coordination number of a grain is large. A 
related scheme, used-to study "grains " with a degenerate 
level is presented intll. This model is a natural extension 
of the Hubbard model used to describe the transition be- 
tween a Mott insulator and a correlated metal, and does 
show this transition. It would be interesting to compare 
the properties of this model to the results presented here. 
A related problem is that of superconducting grains cou- 
pled by the Josephson effect. This model, neglectipg the 
effect of normal quasiparticles, has been treated irc3, us- 
ing the CPA method, which is the analog of the Dynam- 
ical Mean Field Theory used here. 

In the absence of non equilibrium effects, we find that 
the low temperature behavior of the system is insulating. 
There is a well defined Coulomb gap at low temperatures, 
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which, when the intergrain conductances are large, can 
differ significantly from the intragrain charging energy. 
This behavior is different from the low temperature con- 
ductance of a single junction, where the gapless nature 
of the leads gives rise to a power law dependence of the 
conductance, due to cotunneling processes. Our results 
are valid provided that the temperature is much greater 
than the intradot level spacing. Within this restriction, 
and in the absence of non equilibrium effects, we find no 
evidence for a phase transition between an insulating and 
a metallic phasea. 

Non equilibrium effects can lead to a zero temperature 
phase transition, in a similar way to the transition which 
occurs in a single junction. It is interesting to note that 
this phase transition remains restricted to zero tempera- 
ture, even when an array of grains with large coorination 
number is considered. This is so because the most favor- 



able case for a metallic regime is when the environment 
of a given dot is itself metallic. This reduces the prob- 
lem to that of a single grain connected to metallic leads, 
which cannot have a finite temperature phase transition, 
irrespective of the coupling to the leads. 
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FIG. 1: Solution to the large-TV theory for L — 128 and e = 
0.1, 0.3, 0.5, 1.9. 3c decreases with increasing e. Thick 
curves are for e = 0.5 and e = 1.5. 
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FIG. 2: Correlation function G(s) at criticality versus di- 
mensionless imaginary time s. Solid line is input G(s) = 
g-^M + ^4(1 _ e - , 'l s l)| s |-2 e ; dashed line is obtained from Eq. 
|l6| . Here r\ = 3.00 provides the best fit of the three curves. 
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FIG. 4: a c versus e for three different values of ij (curves). 
Open and closed circles denote maximum and minimum val- 
ues of ce c inferred from iterative solution to the self-consistent 
equations with L = 128. 
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FIG. 5: Temperature dependence of the normalized conduc- 
tance g/go, as derived from path-integral Monte Carlo simu- 
lations for a granular system (squares) and for a single tunnel 
junction (circles) with the same parameters: a = 0.1, e = 0, 
and qi = (no offset charge). Error bars are on the order of 
the symbol size. 
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FIG. 6: Normalized conductance (g/go) vs. offset charge q;, 
as derived from Monte Carlo simulations for a granular system 
(squares) and for a single tunnel junction (dashed line) with 
the same parameters: a = 0.1, e = 0, and T — E c /10. 
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FIG. 7: Normalized conductance as a function of a for a 
granular system with e = 0.5. Different symbols represent 
several temperatures: (3E C = 40 (squares), 30 (circles), 20 
(triangles), and 14 (diamonds). 
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FIG. 8: The critical value of the parameter a for the insulator- 
to-conductor transition is plotted vs the parameter e. Sym- 
bols were derived from Monte Carlo simulations, as indicated 
in the text. 



